Spectrum Explorer

In [1]:

    

%matplotlib
import matplotlib.pyplot as plt
import numpy as np
from IPython.html import widgets # Widget definitions
from IPython.display import display, clear_output, HTML # Used to display widgets in the notebook
from IPython.html.widgets import interact, interactive, fixed


    

    




Using matplotlib backend: Qt4Agg

In [2]:

    

from StellarSpectra.model import Model
from StellarSpectra.spectrum import DataSpectrum
from StellarSpectra.grid_tools import SPEX, HDF5Interface
import scipy.sparse as sp
import numpy as np

myDataSpectrum = DataSpectrum.open("../data/Gl51/Gl51RA.hdf5")

myInstrument = SPEX()

myHDF5Interface = HDF5Interface("../libraries/PHOENIX_SPEX_M.hdf5")

myModel = Model(myDataSpectrum, myInstrument, myHDF5Interface, stellar_tuple=("temp", "logg", "Z", "vsini", "vz", "logOmega"), 
                cheb_tuple=("logc0", "c1", "c2", "c3"), cov_tuple=("sigAmp", "logAmp", "l"), region_tuple=("loga", "mu", "sigma"))
myOrderModel = myModel.OrderModels[0]


    

    




Determine Chunk Log: Wl is 2048
Creating OrderModel 0

In [3]:

    

myOrderModel.update_Cov({"sigAmp": 1, "logAmp":-14 , "l":10.})


    

In [4]:

    

#First pass at actually plotting the model
params = {"temp":3000, "logg":5.0, "Z":0.3, "vsini":4, "vz":6, "logOmega":-19.62}
#params = {"temp":4500, "logg":4.0, "Z":-0.1, "vsini":10, "vz":15, "logOmega":-19.}
# params = {"temp":3000, "logg":4.0, "Z":-0.1, "vsini":10, "vz":0, "logOmega":-19.}
myModel.update_Model(params)

#myOrderModel.update_Cheb({"c1":0.0, "c2":0.0, "c3":0.0})
model_flux = myOrderModel.get_spectrum()


    

In [4]:

    

cov = myModel.OrderModels[0].get_Cov()


    

In [5]:

    

spec = myModel.get_data()
wl = spec.wls[0]
fl = spec.fls[0]
mask = spec.masks[0]

# wl_mask = spec.wls[spec.masks]
# fl = spec.fls[spec.masks]

# print(len(wl_mask), len(fl))
# print(wl_mask, fl)
# print(len(wl), len(model_flux))
# print(wl, model_flux)

fig, ax = plt.subplots(nrows=2, figsize=(11,8), sharex=True)
ax[0].plot(wl[mask], fl[mask], "b")
l_model, = ax[0].plot(wl, model_flux, "r")
ax[0].set_ylabel("Data and Model")
l_resid, = ax[1].plot(wl[mask], (fl - model_flux)[mask], "b")
ax[1].set_xlabel(r"$\lambda$\AA")
ax[1].set_ylabel("Residuals")

#cov = myModel.OrderModels[0].get_Cov().todense()

#fig2 = plt.figure()
#ax2 = fig2.add_subplot(111)
#im = ax2.imshow(cov, origin='upper', interpolation='none')

plt.show()


    

In [9]:

    

np.std(fl - myOrderModel.get_spectrum())


    

    
Out[9]:





4.4012350453753493e-16

In [6]:

    

def update_model_plot(**kwargs):
    '''Take the kwargs, update the model and residuals'''
    
    #Update the model spectrum
    myModel.update_Model(kwargs)
    model_flux = myOrderModel.get_spectrum()
    l_model.set_ydata(model_flux)
    
    #Update the residuals
    residuals = (fl - model_flux)[mask]
    l_resid.set_ydata(residuals)
    
    #Find ymax and ymin and rescale
    ax[0].set_ylim(np.min([fl[mask], model_flux[mask]]), np.max([fl[mask], model_flux[mask]]))
    ax[1].set_ylim(np.min((fl - model_flux)[mask]), np.max((fl - model_flux)[mask]))
    
    #Redraw the plot
    fig.canvas.draw_idle()
    
    #Calculate and print the lnprob
    #print(myModel.evaluate())


    

In [11]:

    

np.save("WASPfl.npy", myOrderModel.get_spectrum())
np.save("WASP_resid.npy", fl - model_flux)


    

In [7]:

    

def update_Cheb_plot(**kw):
    '''Take the kwargs, update the model and residuals'''
    
    #Update the Chebyshev polynomial
    myOrderModel.update_Cheb(kw)
    
    model_flux = myOrderModel.get_spectrum()
    l_model.set_ydata(model_flux)
    
    #Update the residuals
    residuals = (fl - model_flux)[mask]
    l_resid.set_ydata(residuals)
    
    #Find ymax and ymin and rescale
    ax[0].set_ylim(np.min([fl[mask], model_flux[mask]]), np.max([fl[mask], model_flux[mask]]))
    ax[1].set_ylim(np.min((fl - model_flux)[mask]), np.max((fl - model_flux)[mask]))
       
    #Redraw the plot
    fig.canvas.draw_idle()
    
    #Calculate and print the lnprob
    #print(myModel.evaluate())


    

In [8]:

    

def update_Cov_plot(**kwargs):
    '''Take the kwargs, update the model and residuals'''
    
    #Update the covariance matrix
    myModel.OrderModels[0].update_Cov(kwargs)
    cov = myModel.OrderModels[0].get_Cov().todense()
    
    
    #Replot the covariance matrix
    im.set_array(cov)
    
    #Redraw the plot
    fig2.canvas.draw_idle()
    
    #Calculate and print the lnprob
    print(myModel.evaluate())


    

In [9]:

    

i = interact(update_model_plot,
         temp=(2900,3200, 10),
         logg=(4,6.0, 0.1),
         Z=(-0.5, 0.5, 0.05),
         #alpha=(0.0, 0.4, 0.05),
         vsini=(3, 8., 0.5),
         vz=(0, 10),
         #Av=(0.0,1.0, 0.05),
         #logOmega=(-19.9,-19.5, 0.01),
         logOmega=(-19.8,-19.3, 0.01),
         )


    

In [13]:

    

i = interact(update_Cheb_plot,
         c1=(-0.2, 0.2, 0.01),
         c2=(-0.2, 0.2, 0.01),
         c3=(-0.2, 0.2, 0.01),
         )


    

    




params are  {'c1': -6.938893903907228e-18, 'c2': -6.938893903907228e-18, 'c3': -6.938893903907228e-18}

In [14]:

    

i = interact(update_Cov_plot,
         sigAmp=(.5,1.5, 0.1),
         logAmp=(-15,-13, 0.2),
         l=(1, 50, 1),
         )


    

    




74024.8847862